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Abstract. Standard assumed displacement finite elements with 
anisotropic material properties perform poorly in complex stress fields such 
as combined bending and shear and combined bending and torsion. To ad- 
dress this problem, a set of three-dimensional hybrid-stress brick elements 
were developed with fully anisotropic material properties. Both eight-node 
and twenty-node bricks were developed based on the symmetry group theory 
of Punch and Atluri. An eight-node brick was also developed using complete 
polynomials and stress basis functions and reducing the order of the result- 
ing stress parameter matrix by applying equilibrium constraints and stress 
compatibility constraints. Here the stress compatibility constraints must be 
formulated assuming anisotropic material properties. The performance of 
these elements was examined in numerical examples covering a broad range 
of stress distributions. The stress predictions show a significant improvement 
over the assumed displacement elements but the calculation time is increased. 
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Introduction 


The development of high strength single crystal metallic alloys, such as those 
used in the turbine blades in the fuel pump in the space shuttle main engine, 
has placed new emphasis on the need for anisotropic stress analysis, especially 
in the area of finite element analysis. These single crystal materials have a 
high degree of anisotropy and the use of standard assiuned displacement finite 
elements can lead to very poor approximations of stresses, displacements, 
natural frequencies, and mode shapes. This work outlines the possibility of 
resolving these deficiencies by developing hybrid stress elements formulated 
for three dimensional linear anisotropic elasticity. 

The hybrid stress elements are based upon a modified complementary 
energy principle in which the displacements and stresses are independently 
interpolated. Two approaches to the interpolations are considered here, both 
of which assure correct stiffness rank, coordinate invariance, and ehmination 
of spurious zero energy modes. The first is based on work by Spilker, Maskeri, 
and Kania [1] and Spilker and Singh [2] in which complete equilibrated poly- 
nomials are used. The number of stress parameters is reduced by applying 
compatibility constraints reformulated in terms of the stresses. The second 
approach is based upon recent work by Rubenstein, Punch and Atluri [3] 
and Punch and Atluri [4, 5] in which group theoretical methods are used 
to minimize the number of stress parameters while still satisfying rank and 
invariance requirements. 

Eight node and 20 node brick elements have been compared to both 
standard displacement elements and to exact solutions, when possible, for 
both isotropic and anisotropic materials. These comparisons have been made 
on single elements as well as on multielement cantilever beams under various 
loading conditions. The hybrid elements have demonstrated significantly 
improved performance for the analysis of highly anisotropic materials. In 
fact, eight node hybrid elements have shown superior performance to 20 
node displacement elements where the material is highly anisotropic. 


1 Derivation of the Element Stiffness Matrix 


The hybrid stress model is developed in the same manner regardless of the 
form of stress interpolation. The development is based upon a modified 
complimentary energy principle in which Lagrange multipliers are used to 
relax interelement traction continuity and mechanical boundary conditions. 
The required functional for the hybrid stress formulation is given by: 

? {5 L - L + L 
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where a is the stress vector, 5 is the compliance matrix for the ^ 

is the strain-displacement matrix, u represents the local displaceme , 

T are the prescribed traction boundary conditions. 

The local displacements are interpolated in terms of the nodal displace- 
ments through the displacement shape functions N. 

u = T7, C)«. 

W = T], C)«^. 

where u. and w, represent the displacements at node i. Equations (1.2) 
can be expressed more concisely in matrix form as 

u = Nq 

where q is a column vector of the nodal displacements. The displacement 
shape functions used for these three dimensional elements are the seren ^p 
itv shape functions. The same shape functions are used to map the ac 
tual element geometry into the master element coordmate system 8en«^^ng 
Standard isoparametric elements in terms of displacements. The coordinate 
transformation is given by 

X — 'Hi 

y = 

^ T?, C)2. 

The strains are related to the displacements in exactly the same manner 
as for a standard assumed displacement technique. The strain displacemen 
matrix D is then given by 


djdx 0 0 

0 djdy 0 
0 0 djdz 

djdy djdx 0 
d/dz d/dz didy 


(1.5) 


and the local strains are given by 


e = Du = DNq = Bq 


( 1 . 6 ) 


The stresses are interpolated separately from the d.sp acements usmg a 
different set of interpolation functions. In addition, rather than interpolating 
"s of the nodJ displacements, the stresses - 

of a set of stress parameters 0 which are only indirectly related to the no 
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displacements. The stresses a are related to these stress parameters via the 
stress interpolation function matrix P: 

a = P(x,y,z)/3 (1.7) 

such that the homogeneous equilibrium conditions are exactly satisfied, i.e., 

Ea = EP^ = 0 (1.8) 


Substituting (1.7) into (1.1) yields 

n™. = E / f^^P^SP(3dV - jT /3'^p'^BqdV + q^N'^TdS^ (1.9) 

The differential volume dV can be found in terms of the determinant of the 
Jacobian of the coordinate transformation as 

dV == \J\didqdQ (1.10) 

The following element matrices can be defined from equations (1.9) and 

(1.10) to simplify statement of the complimentary energy function 

H = P'^SP\J\d^dr]dC 

G = j^j'^j'^P'^B\J\didqdC (1.11) 

Q = [ N'^TdS 

J San 

The element degrees of freedom are related to the global degrees of free- 
dom by the matrix L 

q = Lq- ( 1 . 12 ) 

Substituting equations 1.11-1.16 into equation 1.9 yields: 

Hmc = E (113) 


Taking the first variation of Hmc with respect to /? and setting the varia- 
tion to zero gives: 

{Hl3-GLq*)6(3 = 0 (1.14) 

/3=^H-^GLq^ (1.15) 


Substituting this definition for /? back into yields: 

n.. = E {yirarn-^GL,- - Lror 
= E {-\q*"L^G'^H~^GLq- + q-^LTQ^ 


(1.16) 
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Setting the first variation of lime with respect to q to zero, we obtain 

i-L'^G'^H-^GLq* + L'^Q) 6q* = 0 

^ (1.17) 

L'^KLq* = L'^Q 

where K = G^H~^G and is the element stiffness matrix and Q is the forcing 
vector due to surface tractions. 


2 Formulation of the Assumed Stress Field 

The assumed stress field must satisfy the homogeneous equilibrium equa- 
tions must provide an element stiffness matrix with no spurious zero energy 
modi, and must be invariant under rotation. The number of stress param- 
eters must be greater than or equal to the number of element degrees of 
freedom minus the number of rigid body modes for the stiffness matrix to be 
of sufficient rank to eliminate spurious modes. This minimum number of s 
is a necessary condition but not a sufficient condition to guarantee adequa e 
rank in the stiffness matrix, nor does a stress field which achieves adequate 
rank assure an element which is rotationally invariant. However, systematic 
approaches to defining the stress field are available which can achieve both 

goals. 

Spilker has shown that if the stress field consists of complete polynomials 
(i.e., complete quadratics, or complete cubics) the resulting stiffness matrix 
will both be of sufficient rank and will be invariant so long as the number of 
5’s equals or exceeds the minimum requirement. For an eight node brick a 
full quadratic stress interpolation is required while for the 20 node brick the 
minimum complete polynomial is a cubic. The stress field for the eig no e 
bricks is then defined as; 

Pi 

P2 

P3 

Pa 

Ps 

P6 

where the pi's are row vectors of length 10. Thus 

= ^1+ 2^23: + + 2^rxz -f /3sy^ + 2^9y^ + (2.19) 

A total of 60 /?’s are required to initially define the stress field, but when 
equilibrium constraints are applied, twelve of these stress parameters are 
eliminated. To further reduce the number of stress parameters, the following 
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strain compatibility conditions are applied: 



For a fully anisotropic material, the strains are defined in terms of the com- 
ponents of the compliance matrix as: 


— Cliiax a,12<^y + Ctl3<^z + 0.lATxy + “isT'xz + O-lGTyz (2.21) 

where the <j’s are defined as in (1.18) above. Applying the first compatibility 
constraint yields 

a' , 

2 “1“ ^\2^y d~ ^Xb'^xz d” ^16*^j/z) 


U . 

d” q^2 (^12^37 d“ (^22^ y d“ ^2Z^z d“ U24T^y d“ ^252”xz d“ ^26^yz) (2.22) 


dxdy 


(^14^x d“ (^2A^y d“ ^34^z d” ^44*^xy d” ^45'^xz d” ^462~yz) — 0 


Carrying out the differentiation on the equilibrated stress field yields 
^11 (/?4)d-ai2(~^47d-«13( As d*ai4(/?26) d- 015 (^ 36 ) d-Ctl6(/^46) 
d-Ui 2 (— As) d-a 22 (Aod- 023(^6 -1-024(^3) d-a 25 (/^ 33 )d- 026(^3) 

— (ai4( — A6)d-U24( ~ A 3 ) - 1 - 034 (^ 7 ) - 1 - 044 (^ 4 ) d-a45(/d34)d-a46(A4)) = 0 


(2.23) 


Each equation can be used to eliminate a single /? so that the resulting stress 
field contains 42 stress parameters. 

The same procedure can be applied to the 20 node brick except that 
since the polynomials are cubics, the resulting equations will be functions of 
X, y, and z. Since these compatibility constraints must be independent of 
position, each of these coefficients can be equated to zero yielding a total of 
24 equations. Three of these equations are dependent so that a total of 21 
^’s may be eliminated. Starting with a full cubic stress field with 120 /?’s 
and eliminating 30 stress parameters via the equilibrimn constraints and an 
additional 21 parameters via the compatibility constraints produces a stress 
field with 69 /3's. 

An alternative approach to the stress field is based upon the work of 
Punch and Atluri. The element is formulated as a cube and due to its sym- 
metry, it is invariant to certain rotational and reflective transformations. The 
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, • ^ nf 24 invariant transformations 

symmetry group of a cube comprise rotations six reflections, and 

consisting of one identity trans orma p-j-oun has five irreducible repre- 

eight rotation^reflection. Ttas symmetry h- five 

sentations. The displacement held ^0 node bncK 

quadratics. This displacement field end these 

Xp^es are proj3“lrnaTural irL^^^ 

invariant stress subspaces. 

At this point, the stress field is precisely ‘Je 

stress field produced above. However, smce each ^esult- 

the removal of any stress subspace dojs not idler th ^..^ed from 

ing stress field. An invariant least °*^73«^'thTce" for the least 

54 of ‘h- are invariant and form 
order stress field. Wh th necessarily formed from incomplete 

stiffness matrices of full ran , ^ . , cardinal stress states for pure 

polynomials and thus may no though provide excellent rep- 

oTrarS^ir^s St: es for -sion, ^em 
"‘‘basic premise on which these 

groups for the cube. For amsottop.c ,o the 

violated, but the resulting element is a reason elements but 

stress field. This same lack of symmetry tolerate 

Punch and Atluri have shown that e j/ ovo- distorti^^^ while 
mild distortions and 20 node bricks should have 

still providing good results. The asymmetry due to anisotropy 

similar effects on element performance. 

The elements formed from “7«:rd 

rium and compatibility 7*^^t"th:v°:iso suffer in time comparisons since 
highly anisotropic materials hu y ^ 20 node 

the number of stress parameters is so high (42 and 

bricks respectively). 

3 Determination of the Degree of Anisotropy 

A fully anisotropic material *7 materials properties 

one of these parameters However to determine the influence of 

differ from an isotropic materiah Howeve^^ 

Te::‘s^ftoTy‘s“l3:; Imy the material properties from an isotropic 
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material to a fully anisotropic material. Since very few materials exhibit 
full anisotropy in all orientations, it seems appropriate to begin with an 
orthotropic compliance matrix and by performing a non-symmetric rotation 
on it, generate a fully anisotropic compliance matrix. If the rotation is held 
constant, then by varying the ratio of the Young’s moduli, shear moduli, 
or Poisson’s ratios, compliance matrices varying from fully isotropic to fully 
anisotropic can be generated. For these tests, the following direction cosine 
matrix was used for the rotation: 


.5774 .5774 .5774 ' 

.7071 -.7071 .0000 

- .4082 -.4082 .8165 _ 

Both the ratio of the Young’s moduli and the shear moduli were varied 
over the range of .1 < jE’n/jFaa < 10 and .1 < G 12 /G 23 ^ 10. The remaining 
modulus was set to the average of the other two. While many other schemes 
for systematically varying the degree of anisotropy could be used, this method 
appears adequate for this purpose. 


4 Numerical Results 

The eight node and twenty node hybrid elements were compared to the stan- 
dard displacement elements for both single elements and for a six element 
beam. The single elements were tested in pure tension, pure shear, bending, 
and torsion using both isotropic and anisotropic properties. As is shown in 
Table I, all elements gave exact solutions for pure tension and pure shear 
with isotropic material properties. For pure bending, the eight node dis- 
placement element is overly stiff but all of the hybrid elements again give 
exact results. None of the elements is able to give exact results for torsion, 
but the hybrid elements perform as well or better than the corresponding 
displacement element. 

Table I. Displacements Produced by Cardinal Stress States for 
Isotropic Material Properties 



Pure 

Pure 

Pure 

Pure 

Element 

Tension 

Shear 

Bending 

Torsion 

DM 8 

100 

100 

67 

84 

H8-42 

100 

100 

100 

84 

H8-18 

100 

100 

100 

84 

DM 20 

100 

100 

100 

95 

H20-54 

100 

100 

100 

102 
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As the degree ot anisotropy is increased, the ab.h y of ‘^e f-plac=men^ 
tor anisotropic material properties. Here the ratio c>t 

material axis is rotated with respect to the element axis by the direct^io 
cosines given in the previous section. Both the 

deterioration in bending and torsion the degree of anisotropy mcreas 
though the degradation is small for the 20 node element. 

Table II. Displacements Produced by Cardinal Stress States for 
Anisotropic Material Properties 


Element 

Pure 

Tension 

Pure 

Shear 

Pure 

Bending 

Pure 

Torsion 

DM 8 

100 

100 

46 

76 

H8-42 

100 

100 

100 

84 

H8-18 

100 

100 

100 

84 

DM 20 

100 

100 

97 

92 

H20-54 

100 

100 

100 

95 


These calculations were performed using double precision for all real vari- 
able calculations. The results for the 20 node elements were compared foi 
three integration rules: a 4 X 4 X 4, a 3 X 3 X 3, and the 14 point 
by Irons [6]. The differences in results were 

stated that the 14 point rule produced some ° 

such ill conditioning was detected in these runs. Consequently the U po 
rule was used for all subsequent calculations except 3 

assure that the results were indeed the same for the 14 pom 

rules . 

A six element cantilever beam was analyzed using both eight node and 20 
node bricks. These beams are shown in Figure 1. The beams were analyzed 
boi with a pure moment loading and a uniform end shear and ^r isotropic 
material properties as well as a series of anisotropic material properties. Fig^ 
nre 2 shoL the normalized tip displacement for a pure moment oading wi ^ 
eight node bricks as a function of the degree of anisotropy where the degree 
aidsotropy is given by the ratio of the Young’s moduli m the primary material 
axes Sm f shows the normalized tip displacement of the cantilever beam 
for uniform end shear as a function of the degree of anisotropy. The hybrid 
stress elements are clearly less sensitive to the degree of anisotropy than the 
displacement. This is true even for the least-order formulations of Pune 
and Atluri which depend upon element symmetry for their formulation. 

Figure 4 shows a comparison of a, in the cantilever beam for a pure 
moment load on the end of the beam as a function of the degree of ^^i^otropy 
The stresses in the beam should not change as the material properties chang 
"sses except should be zero. This is cleavly uot the case for 
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the displacement elements, even for the 20-node brick. For cases of severe 
anisotropy, |cTy|,^aa, is as much as 93 percent of \cTy\max while the maximum 
stress of 74 percent of \(Jy\max- For the 20-node hybrid stress element the 
corresponding values are 5.6 percent and 3.2 percent. The eight node element 
actually gives better results than the 20 node brick because the stresses were 
interpolated at the 2x2x2 Gauss points which are the optimum points. 

The hybrid stress elements clearly give better displacements and stresses 
for highly anisotropic material properties than their corresponding displace- 
ment elements but at some calculational expense. The calculation times are 
shown in Table III for calculation of the six element cantilever beam. These 
elements require up to three times as long for the calculations but at least 
twice as many elements are required to obtain the same degree of accuracy 
in and the accuracy in the shear stresses and in cf^ and Gy is still better 
in the hybrid stress elements. 

Table III. Calculation Time for 6 Element Cantilever Beam 

DM 8 DM 20 H8-18 H8-42 H20-54 

Time (sec) 67 552 108 562 1422 

5 Conclusions 

The hybrid stress elements presented here can provide significantly improved 
accuracy in both displacements and stresses for highly anisotropic materials 
in areas of high stress gradients. The 20-node hybrid stress brick element 
provides increased accuracy over the 20-node displacement element at a cost 
of approximately a 3 to 1 increase in computation time. The eight node hy- 
brid element H8:18 provides much improved results over the standard eight 
node displacement element with less than twice the computation time. The 
most surprising result, however, is that the eight node hybrid element pro- 
vides almost the same degree of accuracy as the 20 node displacement element 
at one-fifth of the calculational effort. For high degrees of anisotropy, this 
element gives superior results to the 20-node displacement element. 
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Normalized Tip Displacement 




Q 2 » t 1 1 I 1 ■ i I I I I I,.. — i I I I I I 

.1 1 10 

Young's Modulus Ratio 

Fig. 2 - Normalized Tip Displacement for Cantilever Beam with Pure Moment Loading 
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Normalized Tip Displacement 
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